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Abstract— The task of image registration can be divided into 
two major components, i.e., the extraction of control points 
or features from images, and the search among the extracted 
features for the matching pairs that represent the same feature 
in the images to be matched. Manual extraction of control 
features can be subjective and extremely time consuming, and 
often results in few usable points. On the other hand, automated 
feature extraction allows using invariant target features such as 
edges, corners, and line intersections as relevant landmarks for 
registration purposes. In this paper, we present an extension 
of a recently developed morphological approach for automatic 
extraction of landmark chips and corresponding windows in 
a fully unsupervised manner for the registration of multispec- 
tral images. Once a set of chip-window pairs is obtained, a 
(hierarchical) robust feature matching procedure, based on a 
multiresolution overcomplete wavelet decomposition scheme, is 
used for registration purposes. The proposed method is validated 
on a pair of remotely sensed scenes acquired by the Advanced 
Land Imager (ALI) multispectral instrument and the Hyperion 
hyperspectral instrument aboard NASA’s Earth Observing-1 
satellite. 

I. Introduction 

Many current and future applications in Earth Science, 
Space Science, and Exploration Science (will) require inte- 
gration of data acquired from multiple sources. This type 
of integration provides necessary information for intelligent 
navigation and decision support systems in many real-time 
applications. Depending on the application at hand, the in- 
tegration may be performed on-board or on-the-ground. In 
the case of Earth Science, for example, multiple data sources 
are due to satellite, aircraft, and ground measurement sys- 
tems, where the data represent multiple temporal, spectral, 
and spatial resolutions. For the multisource data integration 
described above, a fast, automatic, reliable, and accurate image 
registration is essential. 

Assuming that image data have been systematically cor- 
rected through prior navigation, image registration (or ”pre- 
cision correction”) corresponds to feature-based matching of 
the image data, which refines georegistration to subpixel 
accuracy. Image registration is thus defined as the process 
that determines the best match between two or more images 
acquired at the same time or later, and wigth the same sensor 
of a different one. One dataset is regarded as the reference 
data and all other data, called input data (or sensed data), are 


matched relative to the reference data. The general process of 
image registration includes three main steps: 

1) Extraction of features to be used in the matching pro- 
cess; 

2) application of a feature matching strategy subject to a 
specific metric; 

3) resampling of the data based on the correspondence 
computed from matched features. 

A large number of automatic image registration methods 
have been proposed and surveys can be found in [1], [2]. In 
particular, previous work has focused on the comparison of 
different choices for Steps 1, 2 above, where either an entire 
scene was registered (against a reference scene), or where 
several (small) chips, extracted typically around characteristic 
image features or landmarks from the reference scene, were 
registered initially against corresponding windows from the 
input scene [3]. The working assumption was that, once a 
registration system became operational, a database of such 
landmark chips would be available for use. Developing the 
capability of extracting these chips automatically and inde- 
pendently of any database presents the following relevant 
advantages: 

• Elimination of a chip database that needs to be maintained 
regularly and adapted to each sensor; 

• extracting automatically chip/windows that contain a very 
small amount of clouds, provided that preprocessing (of 
reference and input data) includes cloud detection; 

• processing images of any size; 

• reducing the initial registration error by extracting simul- 
taneously chips (from the reference image) and corre- 
sponding windows (from the input image), thereby en- 
hancing feature matching (or any optimization technique), 
as far as accuracy and speed are concerned; given that the 
initial conditions are closer to the parameters of the final 
transformation, these techniques serve as ’’refinement 
methods.” 

In this paper, we expand a recently proposed morphological 
method for automatic extraction of reference chips and corre- 
sponding input windows from remotely sensed image data [4]. 
Here, this method is adapted to the specific case of automatic 
registration of multispectral datasets. The paper is organized as 



follows. Section II outlines a previously developed morpholog- 
ical method which performs automated feature extraction for 
registration purposes. Section III addresses the modifications 
introduced in the morphological extraction algorithm in order 
to adapt it to processing of multispectral data sets. Section 
IV provides experimental results to validate the proposed 
approach. Finally, Section V concludes the paper with some 
remarks and hints at plausible future research. 

II. Morphological Feature Extraction 

Mathematical morphology [5] is a nonlinear image process- 
ing technique that was originally established by introducing 
fundamental operators applied to two sets. One set is processed 
by another set with certain spatial properties (i.e., shape and 
size), known as structuring element , which is translated over 
the spatial domain of the image. The structuring element 
acts as a filter for extracting or suppressing specific image 
structures. Following standard notation [5], let us consider a 
grayscale image /, defined on the discrete space Z 2 , and a 
structuring element designed by B € Z 2 . The latter is usually 
‘flat’ in the sense that it is defined in the spatial domain of 
the image (the x-y plane). The morphological erosion of / by 
B is defined as 

( feB)(x,y ) =Min (Sit)eZ 2 (B) /(x + s,t/ + t), (1) 

where (x, y) £ Z 2 and Z 2 (B ) denotes the set of discrete 
spatial coordinates associated with pixels lying within the 
neighborhood defined by B. In contrast, the morphological 
dilation of / by B is defined by 

(/ © B)(x, y) = Max (s , t)eZ 2 (B) /(a; - s,y ~ t). (2) 

Using the same notation above, morphological opening and 
closing filters [5] can be respectively defined by: 

(foB)(x,y) = ((feB)®B)(x,y), (3) 

(f*B){x,y) = ((f®B)eB)(x,y). (4) 

The two operators above have been used successfully to 
build morphological profiles by using combinations of struc- 
turing elements of progressively increased size [6]. However, 
these profiles cannot capture the directional information which 
is crucial to extract features such as edges, comers, and line 
intersections. To address this issue, we resort to so-called 
scale-orientation morphological profiles (SOMPs). 

In order to mathematically define the concept of SOMPs, we 
first denote by B p ^dx,dy) a line segment structuring element 
with p pixels along the line (in this work, we set p = 5) and 
slope dy/dx. Using these notations, we can define the SOMPs 
by opening at a given pixel (x, y) of an image / as 

D(dx,dy)( X ->y) “ {[(/ ° B p ^dx,dy)) e B p ^dx,dy)] (^, 2/)}j (5) 

where p = {1, 2, * • • and k is the maximum number of 
considered scales. In similar terms, we can define the SOMPs 
by closing at /(x, y) as follows: 

D(dx,dy)( X i V ) = {[(/ * Bp,(dx,dy)) 0 B p ,(dx,dy )\{ x 5 ?/)}• (6) 


In both cases, a measure of line strength can be computed 
for each scale and orientation by calculating a point-wise 
distance dist between the pixel /(x,y) in the original image 
and the pixel at the same location in the image filtered by the 
considered line segment structuring element as follows: 


s°{x,y) = {dist(f(x,y),[(foB p , (dx4y) )*B P ' {dX ' dy) }{x,y))}, 

(7) 

s *{x,y) = {dist(f(x,y), [(f»B p ,(dx,dy)) 0 ^p : (dx,dy)] { x i 2/))}* 

( 8 ) 

The resulting values from Eqs. (7) and (8) are combined 
into a feature vector D(x,y) — s °{x,y) U /(x,y) U s* (x, y) 
with dimensionality L = 2k x ra, where k is the number of 
scales and m is the number of orientations. 

In order to automatically extract significant points of interest 
for registration, we use the concept of self-information [7]. 
(The more irregular the SOMP associated with a certain pixel, 
the higher the chance that it represents a comer pixel of an 
object.) Let the L components of D(x, y) be denoted by 
{^(x,y)}^ =1 ; then the self-information provided by the l- 
th component (1 < l < L), can be defined by Ii(x,y) = 
-log pi(x : y ), where 


Pi{z,y) = 


di{x,y) 

Eti d i(x,y) 


(9) 


Using the above definitions, the entropy of D(x,y) can be 
obtained from 


L 

H(D(x,y)) = ~Y^Pi( x >y) x log Pi{x,y). (10) 

1=1 


III. Extension to Multi/Hyperspectral Images 

The SOMP operators described in Section II have been used 
in previous work to automatically extract a set of landmark 
chips from a reference scene, and to find corresponding 
windows in the input scene in a fully unsupervised manner [4J. 
A measure of spectral information divergence (SID) [7] was 
first used to establish a preliminary matching of SOMP-based 
profiles in the input scene with those in the reference scene, 
followed by a (hierarchical) robust feature matching (RFM) 
procedure that makes use of a multiresolution overcomplete 
wavelet decomposition scheme to extract chip and window 
features at the various levels of the decomposition [8]. 

Although the proposed approach was successful for the 
registration of single spectral bands extracted from Landsat 
data sets [4], the following important aspects, which are crucial 
for the registration of multi/hyperspectral image datasets, were 
not addresed: 


• Our previous approach did not incorporate any cloud 
detection procedure and the registration process was 
strongly affected by atmospheric interferes such as 
clouds and their shadows. In order to address this issue, 
we have resorted to spectral unmixing techniques [7] 
to automatically detect clouds (and their shadows) from 
the input scenes. Specifically, we considered the full 
spectral information in the data and extracted the brightest 



pixel vector using the automatic target generation process 
(ATGP) described in [7], [9]. The second iteration of 
ATGP provides the pixel which is most different, in an 
orthogonal sense, from the brightest pixel (essentially, the 
darkest pixel in the scene). These pixels (often associated 
with clouds and (cloud) shadows) are used to mask 
out cloudy/shaded areas by unmixing the scenes using 
linear unmixing [9] and removing pixels with very high 
abundance (90% or more) of either cloud or shadow prior 
to registration. 

• If the pair of images to be registered are rotated against 
each other, as is indeed the case with ALI-Hyperion pairs, 
comparing the SOMPs directly (by using, e.g., a measure 
of spectral similarity) does not make sense, as opposed 
to the registration of Landsat data, for which a trivial 
rotation close to 0 degrees is expected. In fact, the SOMPs 
can be seen as histograms that indicate the degree of 
variation at different scales and orientations (see [6] for 
a more detailed explanation). Thus, if we calculate such 
a histogram for a certain pixel, then rotate the image, 
and recalculate the histogram, the order of the histogram 
bins (which show a measure of change of the pixel with 
regards to surrounding features in the scene at different 
scales and orientations) would be different, even though 
they refer to the same image pixel. To address this issue, 
we simply order the bins in the histogram in descending 
order (i.e., the lowest measures of change are placed first 
and the higher measures of change are placed last in the 
histogram) prior to spectral similarity matching using the 
SID measure. This issue is further discussed in [6]. 

Taking in mind the issues above, the sequence of steps 
implemented by our revised SOMP-based algorithm is: 

1) Remove clouds and shadows in both the reference and 
input multispectral scenes using spectral unmixing. 

2) Obtain the SOMPs for each pixel in the input scene, 
denoted by D g (x,y ), and order the bins in descending 
order for each calculated SOMP. 

3) For each pixel f(x, y) in the reference scene, construct 
its SOMP denoted by Df(x , y) and calculate its entropy 
H(D f (x, y)). 

4) Select the pixel f(x',y f ) with the highest H(Df(x,y)) 
score (i.e., with the maximum entropy) in the reference 
scene and extract a reference chip centered on /(#', y'), 
where the chip size is defined by an input parameter. 

5) Order the bins of the histogram given by Df(x',y') 
in descending order and compute the SID between 
Df(x',y') and all the SOMPs in the input scene ob- 
tained in Step 2, retaining the pixel coordinates (x'\ y") 
of the pixel in the input scene for which the SID value 
is the smallest. 

6) Extract a window centered on pixel (x'\y") in the input 
scene, where the window size is defined by an input 
parameter. 

7) Return to Step 4 until a predefined number of different 
chip-window pairs is extracted. 


TABLE I 

ALI-Hyperion Spectral Band Correspondences. 


ALI 

Hyperion 

Band 

Spectral 

Band 

Average 

number 

range Om) 

number 

wavelength ( nm ) 

1 

0.433 - 0.453 

9 

436.99 

2 

0.45 - 0.515 

14 

487.87 

3 

0.525 - 0.605 

22 

569.27 

4 

0.63 - 0.69 

31 

660.85 

5 

0.775 - 0.805 

44 

793.13 

6 

0.845 - 0.89 

50 

854.18 

7 

1.2- 1.3 

106 

1205.07 

8 

1.55-1.75 

149 

1638.81 

9 

2.08 - 2.35 

207 

2224.03 


TABLE II 

Registration Results Obtained for Five Chip-Window Pairs 
from a Pair of ALI-Hyperion Multi/Hyperspectral Images. 


Chip- window 
pair 

Size 

(pixels) 

Rotation 

(degrees) 

Initial shift 

Initial shift 
f y 

1 

128 x 128 

11.28 

-1.442 

0.153 

2 

96 x 96 

10.92 

1.390 

0.380 

3 

128 x 128 

10.84 

-0.020 

-0.70 

4 

128 x 128 

10.92 

0.370 

-1.11 

5 

128 x 128 

10.69 

-0.720 

0.44 


IV. Experimental Results 

We have tested the fully automated system described in 
section IE using a pair of remotely sensed scenes col- 
lected by the Advanced Land Imager (ALI) multispectral 
instrument and the Hyperion hyperspectral instrument aboard 
NASA’s Earth Observing- 1 satellite. Although the instru- 
ments are mounted on the same platform, their spectral 
resolutions are different. For illustrative purposes, Table I 
shows the correspondences established between the 9 spec- 
tral bands provided by ALI and the corresponding Hype- 
rion bands. A detailed description of band coverage and 
correspondences for ALI-Hyperion is available online (from 
http://eol.usgs.gov/Hyperion_Spectral_Coverage.htm). 

In this study, we have selected a single pair of spectral bands 
from the considered ALI-Hyperion scenes for validation. The 
selected pair, displayed in Fig. 1, consists of ALFs band 7 
and Hyperion’s band 106. Both images are of size 256 x 
3352 pixels. Strong atmospheric interferes (i.e., clouds and 
shadows) can be observed in these ALI and Hyperion images 
depicted in the figure. 

The automatic algorithm described in Section IE was ap- 
plied to the above scenes (with k = 5 scales and m = 8 orien- 
tations) to obtain a set of five chip-window pairs displayed in 
Fig. 2. Due to the narrow width of the ALI-Hyperion scenes 
(only 256 pixels), we had to carefully adjust the desired size 
of the extracted chip-window pairs. For all chip- window pairs, 
the size was set to 128 x 128 pixels, except for pair 2 (in Fig. 2), 
whose size was set to 96 x 96 pixels. In this particular case, the 
proximity to the image border prevented us from extracting a 
128 x 128-pixel window around the image feature provided by 
the morphological operators. As can be observed from Fig. 2, 
all selected chip- window pairs are free from clouds/shadows. 

After applying the hierarchical version of the RFM tech- 




Hyperion window 1 Hyperion window 2 Hyperion window 3 Hyperion window 4 Hyperion window 5 


Fig. 2. ALI chips (first row) and corresponding Hyperion windows (second row) extracted automatically by the proposed SOMP-based method. The size of 
chip-window pairs labeled as 1,3, 4, and 5 is 128 X 128 pixels, while the size of chip-window pair 2 is 96 X 96 pixels. 


nique of [8], we obtained the individual chip-window regis- 
trations given in Table II. From these transformations, we used 
a robust estimation algorithm to derive a global transformation 
of 11.48 degrees ( t x — —160.57 and t y = 152.41), which is 
in good agreement with a previously obtained result using a 
manual procedure, which resulted in a global transformation of 
11.34 degrees (t x — —160.28 and t y = 149.8). This indicates 
that the proposed automatic registration approach can produce 
accurate results without prior knowledge. 

V. Conclusions and Future Research 

In this paper, we have described a fully automated system 
for registration of multispectral image datasets. The system 
is composed of two stages. First, a morphological feature 
extraction technique is used to obtain scale-orientation mor- 
phological profiles (SOMPs) centered on visually significant 
landmarks such as edges, comers, and line intersections. These 
profiles are then used to automatically extract a set of initially 
matched chip-window pairs which are then refined using 
a robust feature matching (RFM) procedure. The combined 
SOMP-RFM technique has been tested in this work using 
a pair of remotely sensed scenes collected by the ALI and 
Hyperion instruments aboard NASA’s EO-1 satellite. Although 
experimental results are very promising, further work is 


required to fully accommodate all the spectral information 
available from the multispectral scenes into the registration. 
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